package ownclasses;

import ij.IJ;

import java.awt.geom.Point2D;

public class Eigensolver2x2 {

	/*
	 * The methods below compute the eigenvalues/eigenvectors for a 2x2
	 * matrix
	 *   | A B |
	 *   | C D |
	 * The implementation was inspired by
	 *   Blinn, Jim: Consider the lowly 2x2 matrix. IEEE Computer Graphics 
	 *   and Applications, 16(2):82-88, 1996.
	 */
	
	public static EigenPair[] realEigenValues2x2(double[][] A) {
		return realEigenValues2x2(A[0][0], A[0][1], A[1][0], A[1][1]);
	}

	public static EigenPair[] realEigenValues2x2(double A, double B, double C, double D) {
		final double R = (A + D) / 2;
		final double S = (A - D) / 2;
		if ((S * S + B * C) < 0) // no real eigenvalues
			return null;
		else {
			double T = Math.sqrt(S * S + B * C);
			double eVal1 = R + T;
			double eVal2 = R - T;
			Point2D eVec1;
			Point2D eVec2;
			if ((A - D) >= 0) {
				eVec1 = new Point2D.Double(S + T, C);
				eVec2 = new Point2D.Double(B, -S - T);
			} else {
				eVec1 = new Point2D.Double(B, -S + T);
				eVec2 = new Point2D.Double(S - T, C);
			}

			EigenPair[] e = new EigenPair[2];
			// put eigenpair with larger eigenvalue up front
			//IJ.log("Eval1: " + eVal1 + " | eval2: " + eVal2);
			if (Math.abs(eVal1) >= Math.abs(eVal2)) {
				e[0] = new EigenPair(eVal1, eVec1);
				e[1] = new EigenPair(eVal2, eVec2);
			} else {
				e[0] = new EigenPair(eVal1, eVec1);
				e[1] = new EigenPair(eVal2, eVec2);
//				throw new Error("switching eigenvalues should never happen!");
			}
			return e;
		}
	}

	/*
	 * EigenPair is a tuple <eigenvalue, eigenvector> and represents
	 * the solution to an eigen problem.
	 */
	
	public static class EigenPair {
		double eival;
		double[] eivec = null;

		public EigenPair(double eival, double[] eivec) {
			this.eival = eival;
			this.eivec = eivec;
		}

		public EigenPair(double eival, Point2D v) {
			this.eival = eival;
			this.eivec = new double[2];
			this.eivec[0] = v.getX();
			this.eivec[1] = v.getY();
		}

		public double getEigenvalue() {
			return this.eival;
		}

		public double[] getEigenvector() {
			return this.eivec;
		}
		
		public String toString() {
			if (eivec == null)
				return "void eigenvalue / eigenvector";
			else {
				String s = "eigenvalue: " + eival + " / eigenvector: (";
				for (double c : eivec) {
					s = s + c + ", ";
				}
				s = s + ")";
				return s;
			}
		}
	}
}
